SUBROUTINE m1cg1 (mvprod,n,x,b,a,xx,bb,two,epsneg,restol,absrel,&
                  iter,imp,iom,imode,mode,w,nw,&
                  bfgsp,pmat0,npmat0,m0,ilm0,nilm0,wlm0,nwlm0,&
                  bfgsb,pmat1,npmat1,m1,ilm1,nilm1,wlm1,nwlm1,&
                  select,izs,rzs,dzs)
  USE constants,                ONLY : dp
  USE io,                       ONLY : ionode

!c
  implicit none
  external mvprod
  logical two
  integer i,j,n,absrel,iter,imp,iom,imode(3),mode,nw,&
          bfgsp,npmat0,m0,nilm0,ilm0(nilm0),nwlm0,&
          bfgsb,npmat1,m1,nilm1,ilm1(nilm1),nwlm1,select,&
          izs(1)
  real rzs(1)
  real(kind=dp) ::  x(n),b(n),xx(n),bb(n),a(1),epsneg,restol,w(nw),&
                    pmat0(n,n),wlm0(nwlm0),&
                    pmat1(n,n),wlm1(nwlm1),dzs(1)
!c
!c----
!c
!c     This routine solves for x the linear system
!c
!c         A x = b,                                               (LS)
!c
!c     by a possibly preconditioned Fletcher-Reeves conjugate gradient
!c     (CG) algorithm. In this system, the matrix A and the vector b are
!c     supposed known. To solve the linear system completely, the matrix
!c     A must be symmetric positive definite. The way A is actually
!c     stored is irrelevant for M1CG1, since the algorithm requires only
!c     products Av of A times various vectors v and since this operation
!c     is supposed to be made in the user-supplied routine MVPROD.
!c
!c     The CG iterations can be preconditioned by a full positive
!c     definite matrix strored in PMAT0 in the standard way (when
!c     bfgsp=1) or by a limited memory BFGS (l-BFGS) matrix stored in
!c     the structure (M0, ILM0, NILM0, WLM0, NWLM0) (when bfgsp=2). The
!c     latter may be particularly appropriate when the dimension N of the
!c     LS is very large. The l-BFGS structure, or even PMAT0, may be
!c     filled in by a previous call to M1CG1. This routine has indeed the
!c     feature (when BFGSB is nonzero) to build a preconditioning matrix
!c     using pairs (Av,v), where the vectors v are the generated
!c     conjugate directions. This preconditioning matrix can then be used
!c     to speed-up another CG run for a system having a matrix not too
!c     different from A.
!c
!c     If two is .true. the routine also solves in parallel the linear
!c     system
!c
!c         A xx = bb,                                            (LSbis)
!c
!c     with the same matrix A and another right hand side bb, using
!c     the conjugate directions generated by the CG iterations for
!c     solving (LS). This feature is useful for the truncated SQP
!c     algorithm in constrained optimization. There is no stopping test
!c     based on the behavior of CG on (LSbis). If (LSbis) is solved
!c     before finding the approximate solution of (LS), xx is modified
!c     by adding terms that should be close to zero.
!c
!c
!c     On entry
!c     ^^^^^^^^
!c
!c       mvprod, external: routine making the products A*v, where v is
!c         an arbitrary vector; it must be defined as follows:
!c           subroutine mvprod (n,a,v,av,izs,rzs,dzs)
!c         where, n (I), a (I), and izs (I), rzs (I), dzs (I) have the
!c         same meanings as below, v(n) (I) is a real(kind=dp) ::  vector
!c         and av(n) (O) is a real(kind=dp) ::  vector containing the
!c         product A*v
!c       n, integer: dimension of the LS to solve
!c       x(n), real(kind=dp) :: : x0, initial guess for x
!c       b(n), real(kind=dp) :: : RHS of (LS)
!c       a(), real(kind=dp) :: : address of the matrix A, it is passed
!c         to MVPROD and is not used in M1CG1
!c       xx(n), real(kind=dp) :: : initial guess for xx (on entry) and
!c         solution of (LSbis) (on return)
!c       bb(n), real(kind=dp) :: : RHS of (LSbis)
!c       two, logical; if .false. only (LS) is solved and xx and bb are
!c         not used; if .true. (LSbis) is also solved by the conjugate
!c         direction algorithm, using the conjugate directions generated
!c         when solving (LS)
!c       epsneg, real(kind=dp) :: : if during CG iterations, a Rayleigh
!c         quotient v'*A*v/v'*v less than a threshold is encountered, the
!c         algorithm stops, since the matrix A is not considered to be
!c         positive definite; epsneg must be nonnegative. The threshold
!c         depends on the value of epsneg and imode(3) (see this
!c         argument).
!c       restol, real(kind=dp) :: : tolerance on the residual; it is used
!c         to decide whether the residual is sufficiently small to
!c         consider that the linear system (LS) is solved by the current
!c         x with the desired precision (see 'absrel' for more
!c         information)
!c       absrel, integer: specifies the meaning of the tolerance 'restol'
!c         = 0: absolute residual test, i.e., if during CG iterations,
!c              the norm of the residual |A*x-b| is less than restol, the
!c              algorithm stops, since (LS) is considered to be solved by
!c              the current x with enough precision
!c         = 1: relative residual test, i.e., if during CG iterations,
!c              the relative norm of the residual |A*x-b|/|A*x0-b| is
!c              less than restol, the algorithm stops, since (LS) is
!c              considered to be solved by the current x with enough
!c              precision
!c       iter, integer: max number of iterations accepted
!c       imp, integer: printing level on the channel io
!c          =0: no printings
!c         >=1: printings on input and output, and error messages
!c         >=2: one line at each CG iteration, giving
!c              - iter: iteration counter
!c              - cost: the cost (1/2) x^T A x - b^T x
!c              - |r| (if absrel=0) or |r|/|r1| (if absrel=1)
!c              - <(r,Pr): angle between the residual r=Ax-b and the
!c                preconditioned residual Pr (will be 0 if bfgsp=0)
!c              - (Av,v)/|v|^2: Rayleigh quotient of A along the
!c                conjugate directions
!c              - conj: cosine of the angle between 2 succesive
!c                residuals (for the scalar product defined by the
!c                preconditioning matrix P)
!c              - alpha: step-size
!c              - |x|: norm of the current approxinate solution x
!c              - S: selection flag; if select >=1, 'S' means that the
!c                pair (Av,v) has been selected, '-' means that it has
!c                not been selected;
!c              if two=.true., "cost" and ("|r|" or "|r|/|r1|") is also
!c              printed for (LSbis)
!c       iom, integer: output channel for printings
!c       imode, integer: input mode of M1CG1
!c         imode(1)=0: the initial iterate x is set to 0
!c                 =1: the initial iterate x is the input variable x
!c         imode(2)=0: the initial iterate xx is set to 0
!c                 =1: the initial iterate xx the input variable xx
!c         imode(3)=0: use epsneg as the threshold of the Rayleigh
!c                     quotient for stopping the CG iterations
!c                 =1: use as threshold of the Rayleigh quotient for
!c                     stopping the CG iterations:
!c                     epsneg*max_Rayleigh_quatient_encoutered
!c       w(nw), real(kind=dp) :: : working area
!c       nw, integer: dimension of w,
!c         must be >= 5*n+m0 if two=.false.
!c         must be >= 6*n+m0 if two=.true.
!c       bfgsp, integer: specifies the type of preconditioning matrix
!c         made available to precondition the current CG run
!c         =0: no preconditioning
!c         =1: a full symmetric positive definite preconditioning matrix
!c             of order n is available in the structure (pmat0,npmat0)
!c         =2: a limited memory BFGS preconditioning matrix is available
!c             in the l-BFGS structure (m0,ilm0,nilm0,wlm0,nwlm0), built
!c             up by a previous called to m1cg1
!c       (pmat0,npmat0): if bfgsp=1, this structure must contain
!c         preconditioning information (not used otherwise)
!c         . pmat0(npmat0), real(kind=dp) :: : must contain an order n
!c           symmetric positive definite matrix approximating the
!c           inverse of A
!c         . npmat0, integer: dimension of pmat0, must be >= n**2
!c       (m0,ilm0,nilm0,wlm0,nwlm0): if bfgsp=2, must contain an l-BFGS
!c         preconditioning structure (not used otherwise)
!c         . m0, integer: number of (y,s) pairs stored in the structure,
!c         . ilm0(nilm0), integer
!c         . nilm0, integer: dimension of ilm0
!c         . wlm0(nwlm0), real(kind=dp) :: : contains the (y,s) pairs
!c         . nwlm0, integer: dimension of wlm0, must be >= m0*(2*n+1)+1
!c       bfgsb, integer: specifies the type of BFGS preconditioning
!c         matrix to build during the CG iterations,
!c         =0: don t build any preconditioning matrix
!c         =1: build a full BFGS preconditioning matrix in the BFGS
!c             structure (pmat1,npmat1); if pmat1(1,1)<=0 the updates
!c             start from scratch, otherwise pmat1 is updated
!c         =2: build a limited memory BFGS preconditioning matrix in the
!c             l-BFGS structure (m1,nys1,jmin1,jmax1,ybar1,sbar1,size1);
!c             if (ilm1(1)<=0) or ((ilm1(1)>0) and (select=1)), the
!c             updates start from scratch
!c       select, in the case when bfgsb=2, it monitors the selection of
!c         the pairs (Av,v) to build the l-BFGS preconditioning matrix,
!c         with the following meanings:
!c         =0: no particular selection, FIFO policy
!c         =1: mexican selection (the pairs are distributed uniformely
!c             according to the iteration counter for a particular CG
!c             run)
!c         =2: by the Rayleigh quotient (it is tried to have pairs with a
!c             Rayleigh quotient distributed as uniformely as possible
!c             on the logarithmic scale)
!c
!c
!c     On return
!c     ^^^^^^^^^
!c
!c       x: if mode = 0, this is the solution found for (LS)
!c       xx: if two=.true. (on entry) and mode = 0 (on return), this is
!c         the solution found for (LSbis)
!c       epsneg: value of the last Rayleigh quotient v'Av/v'v
!c         encountered
!c       restol: if mode=0, restol = |A*x-b| (if absrel=0 on entry) and
!c         restol = |A*x-b|/|A*x0-b| (if absrel=1 on entry)
!c       iter: number of iterations actually performed; iter also serves
!c         as the iteration counter during the run, so that the user can
!c         know the iteration number when MVPROD is called
!c       mode: output mode of M1CG1
!c         =0: normal terminaison (stop on restol)
!c         =1: a dimension argument has a wrong value
!c         =2:
!c         =3: A is probably indefinite (Rayleigh quotient less than
!c             epsneg),
!c         =4: maximum number of iterations reached
!c         =5: 
!c         =6: the initial residual is zero
!c         =7: error in the routine dysave, this may be due to a wrong
!c             value of select, a scalar product (y,s)<0 (rounding
!c             error ?), or a negative value of m1
!c         =8: the large increase of the cost, due to rounding errors
!c       (pmat1,npmat1): if bfgsb=1, will contain the plain BFGS
!c         preconditioning structure (not used otherwise)
!c         . pmat1(npmat1), real(kind=dp) :: : will contain the matrix
!c         . npmat1, integer: dimension of pmat1, must be >= n**2
!c       (m1,ilm1,nilm1,wlm1,nwlm1): if bfgsb=2, will contain the l-BFGS
!c         preconditioning structure (not used otherwise)
!c         . m1, integer: number of (y,s) pairs to store,
!c         . ilm1(nilm1), integer
!c         . nilm1, integer: dimension of ilm1
!c         . wlm1(nwlm1), real(kind=dp) :: : contains the (y,s) pairs
!c         . nwlm1, integer: dimension of wlm1, must be >= m1*(2*n+1)+1
!c
!c
!c     Other arguments
!c     ^^^^^^^^^^^^^^^
!c
!c       izs() [integer], rzs() [real], and dzs() [real(kind=dp) :: ]:
!c         these variables are not used by M1CG1, they are just passed
!c         to the user-supplied routine mvprod
!c
!c-----------------------------------------------------------------------
!c
!c Author: Jean Charles Gilbert, INRIA.
!c
!c Copyright 1998, 1999, 2011, INRIA.
!c
!c M1cg1 is distributed under the terms of the Q Public License version
!c 1.0.
!c
!c This program is distributed in the hope that it will be useful, but
!c WITHOUT ANY WARRANTY; without even the implied warranty of
!c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the Q Public
!c License version 1.0 for more details.
!c
!c You should have received a copy of the Q Public License version 1.0
!c along with this program.  If not, see
!c <http://doc.trolltech.com/3.0/license.html>.
!c
!c-----------------------------------------------------------------------
!c
      integer ii,nq,nv,nr,nrr,nrm,naux,nrho,&
          nybar0,nsbar0,nsize0,&
          nybar1,nsbar1,nsize1,nol1
      real(kind=dp) ::  restol2
!c
!c---- initialization
!c
      mode=0
!c
!c---- initial printings
!c
      if (imp.ge.1) then
          write (iom,911)
          write (iom,'(a/4x,a,i14/4x,a,4x,1pe9.2/4x,a,10x,e9.2)',advance='no') &
            " m1cg1 (Version 1.2, October 2011): entry point",&
            "dimension of the problem (n):", n,&
            "Rayleigh quotient threshold (epsneg):", epsneg,&
            "tolerance on residual (restol):", restol
          if (absrel.eq.0) then
              write (iom,'(a)') " (absolute)"
          else
              write (iom,'(a)') " (relative)"
          end if
          write (iom,'(4x,a,i7/4x,a,i22)') "maximal number of iterations (iter):", iter,&
                                          "printing level (imp):", imp
          if (bfgsp.eq.0) then
              write (iom,901)
          elseif (bfgsp.eq.1) then
              write (iom,902)
          elseif (bfgsp.eq.2) then
              write (iom,903) min(m0,ilm0(1))
          end if
          if (bfgsb.eq.0) then
              write (iom,904)
          elseif (bfgsb.eq.1) then
              write (iom,905,advance='no')
              if (pmat1(1,1).le.0.d0) then
                  write (iom,'(a)') " (from scratch)"
              else
                  write (iom,'(a)') " (from the given matrix pmat1)"
              endif
          elseif (bfgsb.eq.2) then
              write (iom,'(4x,a,i0,a/7x,a)',advance='no') &
                "An l-BFGS preconditioning matrix is built, with ", m1,&
                " (Av,v)-pairs,", "starting from "
              if ((ilm1(1).lt.1) .or. (select.eq.1)) then
                  write (iom,'(a)',advance='no') "scratch"
              else
                  write (iom,'(a)',advance='no') "the matrix on entry"
              end if
              if (select.eq.0) then
                  write (iom,'(a)') ", first-in-first-out selection"
              elseif (select.eq.1) then
                  write (iom,'(a)') ", mexican selection"
              elseif (select.eq.2) then
                  write (iom,'(a)') ", Rayleigh-quotient selection"
              end if
          end if
          if (imode(1).eq.0) then
            write (iom,'(4x,a)') "The starting point is zero"
          else
            write (iom,'(4x,a)') "The starting point is given on entry"
          end if
          if (two) then
              write (iom,'(4x,a)',advance='no') &
                "Two linear systems are solved simultaneously, "
              if (imode(1).eq.0) then
                write (iom,'(6x,a)') "with a zero starting point"
              else
                write (iom,'(6x,a)') &
                  "with a starting point given on entry"
              end if
          end if
          if (imp.ge.2) write (iom,910)
      end if
  901 format (/4x,"Plain CG iterations")
  902 format (/4x,"Preconditionned CG iterations")
  903 format (/4x,"l-BFGS preconditioned CG iterations, formed of ",i0," (Av,v)-pairs")
  904 format (4x,"No preconditioning matrix is built")
  905 format (4x,"A full BFGS preconditioning matrix is built")
  906 format (4x,"An l-BFGS preconditioning matrix is built, with ",i0," (Av,v)-pairs")
  910 format (1x,87("-"))
  911 format (1x,87("*"))
!c
!c---- check the arguments
!c
      if (n.le.0) then
          mode=1
          if (imp.ge.1) then
              write (iom,'(/a,i0,a/)') " >>> m1cg1: inconsistent call, n = ", n," should be > 0"
              write (iom,911)
          end if
          return
      end if
      if (iter.le.0) then
          mode=1
          if (imp.ge.1) then
              write (iom,'(/a,i0,a/)') " >>> m1cg1: inconsistent call, iter = ", iter, " should be > 0"
              write (iom,911)
          end if
          return
      end if
      if (epsneg.lt.0.d+0) then
          mode=1
          if (imp.ge.1) then
              write (iom,'(/a,1pe12.5,a/)') " >>> m1cg1: inconsistent call, epsneg = ", epsneg, " should be >= 0"
              write (iom,911)
          end if
          return
      end if
      if (restol.lt.0.d+0) then
          mode=1
          if (imp.ge.1) then
              write (iom,'(/a,1pe12.5,a/)') " >>> m1cg1: inconsistent call, restol = ", restol, " should be >= 0"
              write (iom,911)
          end if
          return
      end if
      if ((absrel.lt.0) .or. (absrel.gt.1)) then
          mode=1
          if (imp.ge.1) then
              write (iom,'(/a,i0,a/)') " >>> m1cg1: inconsistent call, absrel = ", absrel, " should be in {0,1}"
              write (iom,911)
          end if
          return
      end if
      if (two) then
          if (nw.lt.6*n+m0) then
              mode=1
              if (imp.ge.1) then
                  write (iom,'(/a,i0,a,i0/)') " >>> m1cg1: inconsistent call, nw = ", nw, " should be >= ", 6*n+m0
                  write (iom,911)
              end if
              return
          end if
      else
          if (nw.lt.5*n+m0) then
              mode=1
              if (imp.ge.1) then
                  write (iom,'(/a,i0,a,i0/)') " >>> m1cg1: inconsistent call, nw = ", nw, " should be >= ", 5*n+m0
                  write (iom,911)
              end if
              return
          end if
      end if
!c
!c     --- check imode
!c
      if ((imode(1).lt.0) .or. (imode(1).gt.1)) then
          mode=1
          if (imp.ge.1) then
              write (iom,'(/a,i0,a/)') " >>> m1cg1: inconsistent call, imode(1) = ", imode(1), " should be in {0,1}"
              write (iom,911)
          end if
          return
      end if
      if ((imode(2).lt.0) .or. (imode(2).gt.1)) then
          mode=1
          if (imp.ge.1) then
              write (iom,'(/a,i0,a/)') " >>> m1cg1: inconsistent call, imode(2) = ", imode(2), " should be in {0,1}"
              write (iom,911)
          end if
          return
      end if
      if ((imode(3).lt.0) .or. (imode(3).gt.1)) then
          mode=1
          if (imp.ge.1) then
              write (iom,'(/a,i0,a/)') " >>> m1cg1: inconsistent call, imode(3) = ", imode(3), " should be in {0,1}"
              write (iom,911)
          end if
          return
      end if
!c
!c     --- check the preconditioning matrix, if appropriate
!c
      if (bfgsp.eq.1) then
          if (npmat0.lt.n**2) then
              mode=1
              if (imp.ge.1) then
                  write (iom,'(/a,i0,a,i0/)') " >>> m1cg1: inconsistent call, npmat0 = ", npmat0, " should be >= ", n**2
                  write (iom,911)
              end if
              return
          end if
          do i=1,n
              if (pmat0(i,i).le.0.d0) then
                  mode=1
                  if (imp.ge.1) then
                      write (iom,'(/a,a,i0,a,i0,a/)') " >>> m1cg1: inconsistent call, bfgsp = 1 but ", "pmat0(", i, ",", i, ") is < 0"
                      write (iom,911)
                  end if
                  return
              end if
          end do
      else if (bfgsp.eq.2) then
          if (nilm0.lt.m0+4) then
              mode=1
              if (imp.ge.1) then
                  write (iom,'(/a,i0,a,i0/)') " >>> m1cg1: inconsistent call, nilm0 = ", nilm0, " should be >= ", m0+4
                  write (iom,911)
              end if
              return
          end if
          if (nwlm0.lt.m0*(2*n+1)+1) then
              mode=1
              if (imp.ge.1) then
                  write (iom,'(/a,i0,a,i0/)') " >>> m1cg1: inconsistent call, nwlm0 = ", nwlm0, " should be >= ", m0*(2*n+1)+1
                  write (iom,911)
              end if
              return
          end if
          if (m0.lt.1) then
              mode=1
              if (imp.ge.1) then
                  write (iom,'(/a,a,i0/)') " >>> m1cg1: inconsistent call, bfgsp = 2 but ", "m0 = ", m0
                  write (iom,911)
              end if
              return
          end if
          if (ilm0(1).lt.1) then
              mode=1
              if (imp.ge.1) then
                  write (iom,'(/a,a,i0,a/)') " >>> m1cg1: inconsistent call, bfgsp = 2 but ", "ilm0(1) = ", ilm0(1), " is inappropriate"
                  write (iom,911)
              end if
              return
          end if
!c
!c         --- if the selection procedure was 2, ilm0(5:min(mp,ilm0(1))) are pointers
!c
          j=min(m0,ilm0(1))
          if (ilm0(2).eq.2) then
              ii=j+4
          else
              ii=4
          end if
          do i=3,ii
              if ((ilm0(i).lt.1) .or. (ilm0(i).gt.j)) then
                  mode=1
                  if (imp.ge.1) then
                      write (iom,'(/a,a,i0,a,i0,a/)') &
                        " >>> m1cg1: inconsistent call, bfgsp = 2 but ","ilm0(", i, ") = ", ilm0(i), " is inappropriate"
                      write (iom,911)
                  end if
                  return
              end if
          end do
      end if
!c
!c     --- check the constructed preconditioner, if appropriate
!c
      if (bfgsb.eq.1) then
          if (npmat1.lt.n**2) then
              mode=1
              if (imp.ge.1) then
                  write (iom,'(/a,i0,a,i0/)') " >>> m1cg1: inconsistent call, npmat1 = ", npmat1, " should be >= ", n**2
                  write (iom,911)
              end if
              return
          end if
          if (pmat1(1,1).gt.0.d0) then
              do i=2,n
                  if (pmat1(i,i).le.0.d0) then
                      mode=1
                      if (imp.ge.1) then
                          write (iom,'(/a,a,i0,a,i0,a/)') &
                          " >>> m1cg1: inconsistent call, bfgsb = 1", " and pmat1(1,1) > 0, but pmat1(", i, ",", i, ") is < 0"
                          write (iom,911)
                      end if
                      return
                  end if
              end do
          end if
      else if (bfgsb.eq.2) then
          if (nilm1.lt.m1+4) then
              mode=1
              if (imp.ge.1) then
                  write (iom,'(/a,i0,a,i0/)') " >>> m1cg1: inconsistent call, nilm1 = ", nilm1, " should be >= ", m1+4
                  write (iom,911)
              end if
              return
          end if
          if (m1.lt.1) then
              mode=1
              if (imp.ge.1) then
                  write (iom,'(/a,a,i0/)') " >>> m1cg1: inconsistent call, bfgsb = 2 but ", "m1 = ", m1
                  write (iom,911)
              end if
              return
          end if
          if ((select.lt.0) .or. (select.gt.2)) then
              mode=1
              if (imp.ge.1) then
                  write (iom,'(/a,i0,a/)') " >>> m1cg1: inconsistent call, select = ", select, " should be in {0,1,2}"
                  write (iom,911)
              end if
              return
          end if
!c     write (6,'(a,i0)') "ilm1(1) = ", ilm1(1)
          if (ilm1(1).gt.0) then
              if (ilm1(2).ne.select) then
                  mode=1
                  if (imp.ge.1) then
                      write (iom,'(/a,i0,a,a/)') &
                     " >>> m1cg1: inconsistent call, select = ", select, " differs from the one used to form ", "the l-BFGS matrix to improve"
                      write (iom,911)
                  end if
                  return
              end if
              if (select.gt.0) then
                  j=min(m1,ilm1(1))
!c     write (6,'(a)') "ilm1 = "
!c     write (6,'(i0)') (ilm1(i), i=1,j+4)
                  do i=3,j+4
                      if ((ilm1(i).lt.1) .or. (ilm1(i).gt.j)) then
                          mode=1
                          if (imp.ge.1) then
                              write (iom,'(/a,a,i0,a,i0,a,i0,a/)') &
                   " >>> m1cg1: inconsistent call, bfgsb ", "= 2 and ilm1(1) = ", ilm1(1), " but ilm1(", i, ")=", ilm1(i)," is inappropriate"
                              write (iom,911)
                          end if
                          return
                      end if
                  end do
              end if
          else
              ilm1(2)=select
          end if
!c     write (6,'(a,i0)') "ilm1(2) = ", ilm1(2)
          if (nwlm1.lt.m1*(2*n+1)+1) then
              mode=1
              if (imp.ge.1) then
                  write (iom,'(/a,i0,a,i0/)') " >>> m1cg1: inconsistent call, nwlm1 = ", nwlm1, " should be >= ", m1*(2*n+1)+1
                  write (iom,911)
              end if
              return
          end if
      end if
!c
      nq=1
      nr=nq+n
      nrr=nr
      if (two.eqv..true.) nrr=nr+n
      nv=nrr+n
      nrm=nv+n
      naux=nrm+n
      nrho=naux+n
!c
      nybar0=1
      nsbar0=nybar0+n*m0
      nsize0=nsbar0+n*m0
      nybar1=1
      nsbar1=nybar1+n*m1
      nsize1=nsbar1+n*m1
      nol1=nsize1+1
      restol2=restol*restol
!c
      call m1cga (mvprod,n,x,b,a,xx,bb,two,w(nq),w(nr),w(nrr),w(nv),&
                 w(nrm),epsneg,restol2,absrel,iter,imp,iom,imode,mode,&
                 bfgsp,pmat0,&
                 m0,ilm0(1),ilm0(3),ilm0(4),wlm0(nybar0),&
                 wlm0(nsbar0),wlm0(nsize0),&
                 bfgsb,&
                 pmat1,&
                 m1,ilm1(1),ilm1(3),ilm1(4),ilm1(5),wlm1(nybar1),&
                 wlm1(nsbar1),wlm1(nsize1),wlm1(nol1),select,&
                 w(nrho),w(naux),izs,rzs,dzs)
!c
      restol=sqrt(restol2)
!c
!c     write (6,'(a)') "ilm1 ="
!c     write (6,'(i0)') (ilm1(i),i=1,nilm1)
      return
      end
!c
!c--------0---------0---------0---------0---------0---------0---------0--
!c
      subroutine m1cga (mvprod,n,x,b,a,xx,bb,two,q,r,rr,v,rm,epsneg,&
                        restol2,absrel,iter,imp,iom,imode,mode,&
                        bfgsp,pmat0,&
                        m0,nys0,jmin0,jmax0,ybar0,sbar0,size0,&
                        bfgsb,pmat1,&
                        m1,nys1,jmin1,jmax1,jol,ybar1,sbar1,size1,ol,&
                        select,&
                        rho,aux,izs,rzs,dzs)
!c
      USE constants,    ONLY : dp, pi
      implicit none
      external mvprod
      logical two
      integer n,absrel,iter,imp,iom,imode(3),mode,bfgsb,select,bfgsp,&
              m0,nys0,jmin0,jmax0,&
              m1,nys1,jmin1,jmax1,&
              jol(m1),izs(1)
      real rzs(1)
      real(kind=dp) ::  x(n),b(n),a(1),xx(n),bb(n),q(n),r(n),rr(n),v(n),&
         rm(n),epsneg,restol2,&
         pmat0(n,n),ybar0(n,m0),sbar0(n,m0),size0,&
         pmat1(n,n),ybar1(n,m1),sbar1(n,m1),size1,&
         aux(n),ol(m1),rho(m0),dzs(1)
!c
!c----
!c
!c     This routine solves the linear system (LS) Ax=b, for x,
!c     by the possibly preconditioned Fletcher-Reeves CG algorithm.
!c     The preconditioning can come from l-BFGS updates stored in ybar0,
!c     sbar0 and rho.
!c
!c     The residual is denoted by r = Ax-b and the preconditioning
!c     matrix by P (=+- inverse of A).
!c
!c     The first iterate is supposed to be in x.
!c
!c     Input:
!c
!c       n (integer): dimension of the LS to solve.
!c       a (real(kind=dp) ::  vector): name of the vector storing the
!c           matrix A, when necessary.
!c       iter = max number of iterations
!c       imp: printing levels
!c         <= 0: nothing is printed,
!c          = 1: error, terminal, and stopping messages,
!c          = 2: also, one line per iteration.
!c       bfgsb:
!c         =0: don't build a BFGS preconditioning matrix
!c         =1: build a full BFGS preconditioning matrix in pmat1, if
!c             pmat1(1,1) <= 0., the matrix is re-initialized
!c         =2: build a limited memory BFGS preconditioning matrix
!c       bfgsp:
!c         =0: no BFGS preconditioning
!c         =1: full BFGS preconditioning with pmat0
!c         =2: limited memory BFGS preconditioning
!c
!c     Output:
!c
!c       x (real(kind=dp) ::  vector of dim n): approximate solution of
!c           the (LS).
!c       b (real(kind=dp) ::  vector of dim n): RHS of the (LS).
!c       q (real(kind=dp) ::  vector of dim n): auxiliary vector for
!c           storing Av.
!c       r (real(kind=dp) ::  vector of dim n): updated residual.
!c       v (real(kind=dp) ::  vector of dim n): inner CG directions.
!c       rm (real(kind=dp) ::  vector of dim n): auxiliary vector for
!c           storing Pr.
!c>>>>>>>>>>>> a mon avis inutile, pourrait etre stocke dans q.
!c       epsneg = value of the last Rayleigh quotient v'Av/v'v
!c         encountered (not modified if v=0)
!c       iter = number of iterations
!c
!c----
!c
!c --- parameter
!c
      !real(kind=dp) ::  pi
      !parameter (pi=3.1415927d0)
!c
!c --- local variables
!c     restart=.true.:
!c           restarts are allowed. A restart of the CG iterations
!c           occurs at iterations n, 3n/2, 9n/4, ... or when
!c           (Pr_,r)>=0.2|Pr_||r| (Powell's criterion).
!c     r2i = |initial residual|^2
!c     v2 = |v|^2
!c     pprr = previous prr
!c     prr = <Pr,r>
!c     r2 = |r|^2
!c     avv = <Av,v>
!c     rel_r=(|r|/|r1|)**2: relative r
!c     cost = (Ax,x)/2 - (b,x)
!c     x2 = |x|^2
!c
      logical restart
      integer miter,iterr,niter,biter,i,j,ysmode
      character sc
      real(kind=dp) ::  r2i,rr2i,v2,pprr,a1,prr,r2,rr2,avv,rel_r,alpha,&
         aalpha,beta,cost,ccost,mcost,rayl,raylmax,gperp,rm2,dangle,x2,&
         t,sp,nu
!c
!c --- initialisation
!c     miter = max iterations
!c     niter = iteration nbr of the next forced restart
!c     iterr = iteration nbr of the previous restart
!c     biter = bad iteration nbr (will produce failure)
!c
      miter=iter
      niter=n
      iterr=1
      biter=int(sqrt(dble(n*n*n)))
      iter=0
      rel_r=1.d0
      gperp=0.d0
      sc=' '
      if (select.eq.0) sc='S'
      nu=epsneg
      raylmax=0.d0
!c
!c     --- In my experience, it is better not doing restarts
!c
      restart=.false.
!c
!c     --- force initialization of the l-BFGS matrix if appropriate
!c
!c     nys1=0
!c     write (6,'(a,i0)') "nys1 = ", nys1
!c
!c     --- initial residual r = Ax-b (and rr = A xx - bb, if appropriate)
!c         r2i = |initial residual|^2
!c
      if (imode(1).eq.0) then
          do i=1,n
              x(i)=0.d0
              a1=-b(i)
              r(i)=a1
              rm(i)=a1
          enddo
      else
          call mvprod (n,a,x,r,izs,rzs,dzs)
          do i=1,n
              a1=r(i)-b(i)
              r(i)=a1
              rm(i)=a1
          enddo
      end if
      if (two) then
          if (imode(2).eq.0) then
              do i=1,n
                  xx(i)=0.d0
                  rr(i)=-bb(i)
              enddo
          else
              call mvprod (n,a,xx,rr,izs,rzs,dzs)
              do i=1,n
                  a1=rr(i)-bb(i)
                  rr(i)=a1
              enddo
          end if
      end if
!c
      call deuclid (n,r,r,r2i,izs,rzs,dzs)
      r2=r2i
      if (r2i.eq.0.d0) then
          mode=6
          if (imp.gt.0) write (iom,900)
          goto 1000
      end if
  900 format(/4x,"m1cga-WARNING: zero initial residual")
      if (two) then
          call deuclid (n,rr,rr,rr2i,izs,rzs,dzs)
          if (rr2i.eq.0.d0) rr2i=1.d0
      end if
!c
!c     --- rm = Pr
!c         prr = <Pr,r>
!c
!c>>>>>>>>> les 3 lignes ci-dessous me semblent inutiles (4/3/99)
      do i=1,n
          rm(i)=r(i)
      enddo
      if (bfgsp.eq.1) then
          do i=1,n
              rm(i)=0.d0
          enddo
          do j=1,n
              a1=r(j)
              do i=1,n
                  rm(i)=rm(i)+pmat0(i,j)*a1
              enddo
          enddo
      elseif (bfgsp.eq.2) then
          call dbfgsl (n,rm,m0,nys0,jmin0,jmax0,ybar0,sbar0,rho,size0,izs,rzs,dzs)
      end if
      call deuclid (n,rm,rm,rm2,izs,rzs,dzs)
      call deuclid (n,r,rm,prr,izs,rzs,dzs)
!c
!c     --- initial search direction v
!c         v2 = |v|^2
!c
      do i=1,n
          v(i)=-rm(i)
      enddo
      v2=rm2
!c
!c --- cost
!c     cost = current cost
!c     mcost = min cost encountered
!c
      cost=0.d0
      ccost=0.d0
      mcost=0.d0
!c
!c --- printings
!c
      if (imp.ge.2) then
          write (iom,901) sqrt(r2i)
          if (two) write (iom,902) sqrt(rr2i)
      end if
  901 format (4x,"m1cga: initial residual |r1| =",1pe13.6)
  902 format (24x,"LS2 |r1| =",1pe13.6)
!c
!c --- start the CG iterations at 100, end at 1000
!c
  100 iter=iter+1
!c
!c     --- check iter (force convergence, however)
!c
      if (iter.gt.miter) then
          mode=4
          iter=iter-1
          if (imp.gt.0) write (iom,910)
          goto 1000
      end if
  910 format(/4x,"m1cga: max number iterations reached")
!c
!c     --- too many iterations ?
!c
!c     if (iter.gt.biter) then
!c         mode=5
!c         if (imp.gt.0) write (iom,911) biter
!c         goto 1000
!c     end if
!c 911 format(/4x,"m1cga-ERROR: number of iterations exceeds",i6)
!c
!c     --- restart ?
!c
!c>>>>>>>>>> this part of the code has not been modified to take care of
!c           the case when there are two linear systems to solve
!c           simultaneously (two=.true.)
!c
      if (restart.and.(gperp.gt.0.2d0.or.iter.gt.niter)) then
          if (iter.gt.niter) niter=niter*3/2
          iterr=iter
          if (imp.ge.2) write (iom,912)
  912     format(4x,"m1cga: restart")
!c
!c         --- recompute the residuals r and rm
!c
          call mvprod (n,a,x,r,izs,rzs,dzs)
          do i=1,n
              a1=r(i)-b(i)
              r(i)=a1
              rm(i)=a1
          enddo
          call deuclid (n,r,r,r2,izs,rzs,dzs)
          if (bfgsp.eq.1) then
              do i=1,n
                  rm(i)=0.d0
              enddo
              do j=1,n
                  a1=r(j)
                  do i=1,n
                      rm(i)=rm(i)+pmat0(i,j)*a1
                  enddo
              enddo
          elseif (bfgsp.eq.2) then
              call dbfgsl (n,rm,m0,nys0,jmin0,jmax0,ybar0,sbar0,rho,size0,izs,rzs,dzs)
          end if
          call deuclid (n,rm,rm,rm2,izs,rzs,dzs)
          call deuclid (n,r,rm,prr,izs,rzs,dzs)
!c
!c         --- and the search direction v
!c
          do i=1,n
              v(i)=-rm(i)
          enddo
          v2=rm2
      end if
!c
!c     --- q = Av
!c         avv = <v,q> = <Av,v>
!c
      call mvprod (n,a,v,q,izs,rzs,dzs)
      call deuclid (n,v,q,avv,izs,rzs,dzs)
!c
!c     --- check positive definiteness
!c
      if (v2.ne.0.d0) then
          rayl=avv/v2
          raylmax=max(raylmax,rayl)
          if (imode(3).eq.1) nu=epsneg*raylmax
      end if
      if (avv.le.nu*v2) then
          mode=3
          if (v2.ne.0.d0) epsneg=avv/v2
          if (imp.gt.0.and.v2.ne.0.d0) write (iom,913) rayl
          goto 1000
      end if
  913 format(/4x,"m1cga-ERROR: non positive definite matrix,"," <Av,v>/|v|^2 = ",1pe9.2)
!c
!c     --- update the preconditioning matrix when appropriate
!c         y=q=Av, s=v, ys=avv
!c
      if (bfgsb.eq.1) then
!c
!c         --- full inverse BFGS update
!c
!c             --- initialize the matrix pmat1 ?
!c
          if (iter.eq.1.and.pmat1(1,1).le.0.d0) then
!c             call deuclid (n,q,q,a1,izs,rzs,dzs)
!c             a1=a1/avv
              call deuclid (n,v,v,a1,izs,rzs,dzs)
              a1=a1/avv
              if (imp.ge.2) write (iom,914) a1
              do i=1,n
                  do j=1,n
                      pmat1(i,j)=0.d0
                  enddo
                  pmat1(i,i)=a1
              enddo
          end if
!c
!c             --- update pmat1
!c
          call dbfgsi (n,pmat1,q,v,avv,aux,izs,rzs,dzs)
      elseif (bfgsb.eq.2) then
!c
!c         --- limited memory inverse BFGS update
!c
          sc='-'
          call dysave (n,q,v,avv,m1,nys1,jmin1,jmax1,ybar1,sbar1,select,iter,ol,jol,nu,size1,ysmode,4,iom,izs,rzs,dzs)
          if (ysmode.eq.0) sc='S'
          if (ysmode.gt.0) then
              mode=7
              goto 1000
          end if
      end if
  914 format(4x,"m1cga: initialization of the BFGS matrix with OL"," factor =",1pe13.6)
!c
!c     --- new iterate x and residual r
!c         new iterate xx and residual rr
!c
      alpha=prr/avv
      x2=0.d0
      do i=1,n
          t=x(i)+alpha*v(i)
          x(i)=t
          x2=x2+t*t
          r(i)=r(i)+alpha*q(i)
      enddo
      if (two) then
          call deuclid (n,rr,v,t,izs,rzs,dzs)
          aalpha=-t/avv
          do i=1,n
              xx(i)=xx(i)+aalpha*v(i)
              rr(i)=rr(i)+aalpha*q(i)
          enddo
      end if
!c
!c --- stop if the cost increases
!c
      mcost=min(cost,mcost)
      cost=0.d0
      do i=1,n
          cost=cost+(r(i)-b(i))*x(i)
      enddo
      cost=0.5d0*cost
      if (cost.gt.0.999999d0*mcost) then
          mode=8
          if (imp.gt.0) write (iom,915)
          goto 1000
      end if
  915 format(/4x,"m1cga-ERROR: cost increase (rounding error)")
!c
      if (two) then
          ccost=0.d0
          do i=1,n
              ccost=ccost+(rr(i)-bb(i))*xx(i)
          enddo
          ccost=0.5d0*ccost
      end if
!c
!c     --- r2=|r|^2
!c
      call deuclid (n,r,r,r2,izs,rzs,dzs)
      if (two) call deuclid (n,rr,rr,rr2,izs,rzs,dzs)
!c
!c     --- prepare to check conjugacy
!c         gperp = (rmo,r)/|rmo|/|r| == 0 ?
!c
      call deuclid (n,rm,r,gperp,izs,rzs,dzs)
      t=sqrt(r2*rm2)
      if (t.ne.0.d0) then
          gperp=gperp/t
      else
          gperp=0.d0
      end if
!c
!c     --- new rm
!c         rm2 = |rm|^2
!c         new pprr=prr
!c         prr=<Pr,r>
!c
      do i=1,n
          rm(i)=r(i)
      enddo
      if (bfgsp.eq.1) then
          do i=1,n
              rm(i)=0.d0
          enddo
          do j=1,n
              a1=r(j)
              do i=1,n
                  rm(i)=rm(i)+pmat0(i,j)*a1
              enddo
          enddo
      elseif (bfgsp.eq.2) then
          call dbfgsl (n,rm,m0,nys0,jmin0,jmax0,ybar0,sbar0,rho,size0,izs,rzs,dzs)
      end if
      call deuclid (n,rm,rm,rm2,izs,rzs,dzs)
      pprr=prr 
      call deuclid (n,r,rm,prr,izs,rzs,dzs)
!c
!c     --- some printings
!c
      if (imp.ge.2) then
          call deuclid (n,rm,rm,sp,izs,rzs,dzs)
          dangle=0.d0
          if (sqrt(sp*r2).ne.0.d0) dangle=acos(prr/sqrt(sp*r2))/pi*180.d0
          if (mod(iter,50).eq.1) then
              write (iom,'(/4x,a,6x,a)',advance='no') "m1cga: iter", "cost"
              if (absrel.eq.0) then
                  write (iom,'(9x,a,5x)',advance='no') "|r|"
              else
                  write (iom,'(7x,a,2x)',advance='no') "|r|/|r1|"
              end if
              write (iom,'(a,4x,a,3x,a,4x,a)',advance='no') "<(r,Pr)(Av,v)/|v|^2", "conj", "alpha", "|x|"
              if (bfgsb.eq.2) then
                  write (iom,'(3x,a)') "S"
              else
                  write (iom,'()')
              end if
          end if
          write (iom,'(4x,"m1cga:",i5,1x,1pe13.6)',advance='no') iter, cost
          if (absrel.eq.0) then
              write (iom,'(1x,1pe13.6)',advance='no') sqrt(r2)
          else
              write (iom,'(1x,1pe13.6)',advance='no') sqrt(r2/r2i)
          end if
          write (iom,'(1x,0pf5.2,1x,1pe13.6,1x,e9.2,1x,e7.0,1x,e8.1)',advance='no') & 
          dangle, rayl, gperp, alpha, sqrt(x2)
          if (bfgsb.eq.2) then
              write (iom,'(1x,a1)') sc
          else
              write (iom,'()')
          end if
!c         write (iom,917) iter,cost,sqrt(r2/r2i),dangle,rayl,gperp,
!c    &                   alpha,sqrt(x2),sc
          if (two) write (iom,918) ccost,sqrt(rr2/rr2i)
!c         call flush (io)
      end if
!c 916 format (/4x,"m1cga: iter",6x,"cost",7x,"|r|/|r1|",2x,"<(r,Pr)",
!c    &         "(Av,v)/|v|^2",4x,"conj",3x,"alpha",4x,"|x|",3x,"S")
!c 917 format (4x,"m1cga:",i5,1x,1pe13.6,1x,e13.6,1x,0pf5.2,1x,1pe13.6,
!c    &        1x,e9.2,1x,e7.0,1x,e8.1,1x,a1)
  918 format (6x,"LS2 ",6x,1pe13.6,1x,e13.6)
!c
!c     --- stop on the relative residual (restol2)
!c
      if (absrel.eq.0) then
          if (r2.le.restol2) then
              mode=0
              if (imp.gt.0) write (iom,'(/4x,a,1pe9.2,a,e9.2)') "m1cga: stopping criterion satisfied, |r| = ", sqrt(r2)," <= ", sqrt(restol2)
              goto 1000
          end if
      else
          rel_r=r2/r2i
          if (rel_r.le.restol2) then
              mode=0
              if (imp.gt.0) write (iom,919) sqrt(rel_r),sqrt(restol2)
              goto 1000
          end if
  919     format(/4x,"m1cga: stopping criterion satisfied, |r|/|r1| = ",1pe9.2," <= ",d9.2)
      end if
!c
!c     --- FR beta, new v, v2=|v|^2
!c
   80 beta=prr/pprr
      do i=1,n
          v(i)=-rm(i)+beta*v(i)
      enddo
      call deuclid (n,v,v,v2,izs,rzs,dzs)
!c
!c     --- loop
!c
      goto 100
!c
!c --- terminaison
!c
 1000 continue
!c     if (imp.eq.1) write (iom,920) iter,sqrt(r2i),sqrt(r2)
!c 920 format(4x,"m1cga: niter=",i5,", |r|:",d13.6," -> ",d13.6)
!c
!c1001 continue
      if (absrel.eq.0) then
          restol2=r2
      else
          restol2=r2/r2i
      end if
      if (imp.gt.0) then
          write (iom,921)
          write (iom,'(a,i8,1x)',advance='no') " m1cg1: output mode is ", mode
          select case (mode)
              case (0)
                  write (iom,'(a)') "(stopping criterion satisfied)"
              case (1)
                  write (iom,'(a)') "(a dimension argument has a wrong value)"
              case (3)
                  write (iom,'(a)') "(the matrix migth be indefinite)"
              case (4)
                  write (iom,'(a)') "(max iterations reached)"
              case (6)
                  write (iom,'(a)') "(the starting point is a solution)"
              case (7)
                  write (iom,'(a)') "(error in safeguading (Av,v) pairs)"
              case (8)
                  write (iom,'(a)') "(the cost increases too much)"
          end select
          write (iom, '(4x,a,i4/4x,a,2x,1pe22.15/4x,a,2x,1pe22.15/4x,a,2x,1pe22.15)') &
           "number of iterations = ", iter,&
           "initial residual     = ", sqrt(r2i),&
           "final residual       = ", sqrt(r2),&
           "relative residual    = ", sqrt(r2/r2i)
          if (imp.gt.1) then
              call mvprod (n,a,x,q,izs,rzs,dzs)
              cost=0.d0
              do i=1,n
                  cost=cost+x(i)*(q(i)/2.d0-b(i))
!c                 q(i)=q(i)-b(i)
              enddo
!c             call deuclid (n,q,q,t,izs,rzs,dzs)
!c             t=dsqrt(t)
!c             write (iom,923) iter,cost,sqrt(restol2),t
              write (iom,'(4x,a,10x,1pe22.15)') "cost function =", cost
!c
              if (two) then
                  call mvprod (n,a,xx,q,izs,rzs,dzs)
                  cost=0.d0
                  do i=1,n
                      cost=cost+xx(i)*(q(i)/2.d0-bb(i))
                      q(i)=q(i)-bb(i)
                  enddo
                  call deuclid (n,q,q,t,izs,rzs,dzs)
                  write (iom,924) cost,sqrt(t/rr2i),dsqrt(t)
              end if
          end if
          write (iom,925)
      end if
  921 format (1x,87("-"))
  922 format (4x,"number of iterations = ",i4/4x,"relative residual = ",5x,1pe22.15)
  923 format (4x,"number of iterations = ",i4/4x,"cost function =",10x,1pe22.15/4x,"relative residual = ",5x,d22.15/4x,"residual = ",14x,d22.15)
  924 format (/" m1cg1: LS2"/4x,"cost function =",10x,1pe22.15/4x,"relative residual = ",5x,d22.15/4x,"residual = ",14x,d22.15)
  925 format (1x,87("*"))
      return
      end
!c
!c--------0---------0---------0---------0---------0---------0---------0--
!c
      subroutine dbfgsi (n,mat,y,s,ys,smmy,izs,rzs,dzs)
!c
      USE constants,    ONLY : dp
      integer n,izs(1)
      real rzs(1)
      real(kind=dp) ::  mat(n,n),y(n),s(n),ys,smmy(n),dzs(1)
!c
!c----
!c
!c     This routine updates the matrix MAT by the inverse BFGS formula.
!c
!c     Input:
!c
!c       n: order of the matrix
!c       mat: original matrix
!c       y: vector y
!c       s: vector s
!c       ys = (y,s)
!c
!c     Ouput:
!c
!c       mat: updated matrix
!c
!c     Working zone:
!c
!c       smmy(n)
!c
!c----
!c
!c --- local variables
!c
      integer i,j
      real(kind=dp) ::  r,d
!c
!c --- compute (s-My)/(y,s)
!c
      do i=1,n
          smmy(i)=s(i)
      enddo
      do j=1,n
          r=y(j)
          do i=1,n
              smmy(i)=smmy(i)-mat(i,j)*r
          enddo
      enddo
      do i=1,n
          smmy(i)=smmy(i)/ys
      enddo
!c
!c --- r=(s-My,y)/(y,s)^2
!c
      call deuclid (n,y,smmy,d,izs,rzs,dzs)
      r=d/ys
!c
!c --- update M
!c
      do i=1,n
          do j=1,n
              mat(i,j)=mat(i,j)+smmy(i)*s(j)+s(i)*smmy(j)-r*s(i)*s(j)
          enddo
      enddo
!c
      return
      end
!c
!c--------0---------0---------0---------0---------0---------0---------0--
!c
      subroutine dbfgsl (n,d,m,nys,jmin,jmax,ybar,sbar,rho,size,izs,rzs,dzs)
!c
      USE constants,    ONLY : dp
      integer n,m,nys,jmin,jmax,izs(1)
      real rzs(1)
      real(kind=dp) ::  d(n),ybar(n,m),sbar(n,m),rho(m),size,dzs(1)
!c
!c----
!c
!c     compute the product H u, where
!c     . H is the matrix that could be obtained from the m pairs
!c       (ybar,sbar) using the inverse BFGS formula from the diagonal
!c       matrix (size) * (Identity matrix),
!c     . u is a vector of dimension n.
!c
!c     The vector d contains
!c
!c       u (on entry), H u (on output).
!c
!c     rho(m) is a working zone.
!c
!c     Input:
!c
!c----
!c
!c --- local variables
!c
      integer i,j,jfin,jp
      real(kind=dp) ::  r
!c
!c --- return if there is no pair (y,s) in ybar and sbar
!c
      if (nys.eq.0) return
!c
!c --- set jfin
!c
      jfin=jmax
      if (jfin.lt.jmin) jfin=jmax+m
!c
!c --- backward sweep
!c
      do j=jfin,jmin,-1
          jp=j
          if (jp.gt.m) jp=jp-m
          call deuclid (n,d,sbar(1,jp),r,izs,rzs,dzs)
          rho(jp)=r
          do i=1,n
              d(i)=d(i)-r*ybar(i,jp)
          enddo
      enddo
!c
!c --  preconditioning
!c
      r=size
      do i=1,n
          d(i)=d(i)*r
      enddo
!c
!c --- forward sweep
!c
      do j=jmin,jfin
          jp=j
          if (jp.gt.m) jp=jp-m
          call deuclid (n,d,ybar(1,jp),r,izs,rzs,dzs)
          r=rho(jp)-r
          do i=1,n
              d(i)=d(i)+r*sbar(i,jp)
          enddo
      enddo
      return
      end
!c
!c--------0---------0---------0---------0---------0---------0---------0--
!c
      subroutine dysave (n,y,s,ys,m,nys,jmin,jmax,ybar,sbar,select,&
                         iiter,ol,jol,nu,size,mode,plev,iom,izs,&
                         rzs,dzs)
!c
      USE constants,    ONLY : dp
      integer n,m,select,nys,jmin,jmax,iiter,jol(m),mode,plev,iom,izs(1)
      real rzs(1)
      real(kind=dp) ::  y(n),s(n),ys,ybar(n,m),sbar(n,m),ol(m),nu,size,dzs(1)
!c
!c----
!c
!c     Save a pair (y,s)/sqrt(ys) in YBAR and SBAR and update the
!c     Oren-Luenberger factor OLFACT.
!c
!c     Input:
!c
!c       ys: Euclidean scalar product (y,s). Must be positive.
!c       nys: number of (y,s) pairs having been received. Must be
!c         initialized to 0 before calling dysave for the first time.
!c       select monitors the selection of the pairs (y,s) to build the
!c         l-BFGS preconditioning matrix, with the following meanings:
!c         =0: no particular selection, FIFO policy
!c         =1: mexican selection (the pairs are distributed uniformely
!c             on the iteration counter iiter for a particular CG run,
!c             at each new run the previous pairs are discarded and a
!c             new l-bfgs matrix is build from scratch)
!c         =2: selection by the Rayleigh quotient; the aim is to obtain
!c             uniformely distributed OL factors in the logarithmic
!c             scale
!c       iiter: integer giving the index of the current inner iteration
!c         when dysave is called.
!c       jol(m):
!c         . if select = 1, jol(i) gives the iteration number of the i-th
!c           selected pair
!c         . if select = 2, jol(i) is the index n ol of the pair (y,s)
!c           with the i-th smallest OL factor.
!c       nu = epsilon for detecting negative curvature directions,
!c         whose square root is used for selecting good pairs.
!c       size: scalar preconditioner.
!c
!c     Output:
!c
!c       mode: output mode
!c         =-1: the (y,s) pair is not selected
!c         =0 : good call,
!c         =1 : the number m of pairs that can be stored in (ybar,sbar)
!c              is <= 0.
!c         =2 : ys <= 0.
!c         =3 : unknown value of select
!c
!c----
!c
!c --- local variables
!c
      integer i,i2,j,jcour,odiff
      real(kind=dp) ::  rmin,r,olf,olmed,oldq,newq,ol0,ol1,ol2,ol3,olx,olvar,olvar1,olvar2
!c
!c --- check the input
!c
      if (m.le.0) then
          mode=1
          if (plev.ge.1) write (iom,900) m
          return
      end if
  900 format (/" dysave-ERROR: non positive number of updates m = ",i5)
!c
      if (ys.le.0.d0) then
          mode=2
          if (plev.ge.1) write (iom,901) ys
  901     format (/" dysave-ERROR: non positive scalar product"," (y,s) = ",1pe12.5)
          return
      end if
!c
!c --- initialization
!c
      mode=0
      rmin=1.d-20
!c
!c --- Oren-Luenberger factor
!c
!c     call deuclid (n,y,y,olf,izs,rzs,dzs)
!c     olf=ys/olf
!c
!c     for a quadratic function, I prefer
!c
      call deuclid (n,s,s,olf,izs,rzs,dzs)
      olf=olf/ys
!c>          print *,"olf=",olf
      if (plev.ge.5) write (iom,902) olf
  902 format (/4x,"dysave: inverse Oren-Luenberger factor: ","|s|^2/(y,s) = ",1pe10.3)
!c
!c --- select the pair to discard if any and update the pointers
!c
!c     --- is this a good pair to save ?
!c
      if (select.ge.0 .and. olf.lt.sqrt(nu)) then
          mode=-1
          return
      end if
!c
!c     --- accept any pair, discard the oldest one
!c
      if (select.eq.0) then
          if (nys.le.0) then
              nys=0
              jmin=1
              jmax=0
          end if
          nys=nys+1
          jmax=jmax+1
          if (jmax.gt.m) jmax=jmax-m
          if (nys.gt.m) then
              jmin=jmin+1
              if (jmin.gt.m) jmin=jmin-m
          end if
          jcour=jmax
          ol(jcour)=olf
!c         write (6,'(a,i0,a,i0,a,i0,a,i0)') "nys=", nys, "  jmin=",
!c    &      jmin, " jmax=", jmax, "  jcour=", jcour
!c
!c     --- mexican selection
!c
      elseif (select.eq.1) then
          if (iiter.eq.1) nys=0
          if (nys.le.0) then
              nys=0
              jmin=1
              jmax=0
              do i=1,m
                  jol(i)=0
              enddo
          end if
          nys=nys+1
          if (nys.le.m) then
              jmax=jmax+1
              jol(jmax)=iiter
              jcour=jmax
          elseif (m.eq.2) then
              jol(2)=iiter
              jcour=2
          elseif (m.gt.2) then
!c
!c         --- reject current pair if iiter is too close to jol(m-1)
!c             jol gives the iteration number of the selected pairs
!c             they are supposed to be in order of iteration
!c             here jmin=1 and jmax=m (always)
!c
              odiff=jol(m)-jol(m-1)
              if (iiter-jol(m) .lt. odiff) then
                  mode=-1
                  return
              end if
!c
!c         --- select the pair jcour to discard
!c
              do j=1,m-1
                  jcour=j+1
                  if (jol(jcour)-jol(j) .eq. odiff) goto 10
              enddo
   10         continue
!c
!c         --- shift the other (ybar,sbar) and jol
!c>>>>>>>>>>>> playing with pointers instead of with the actual pairs
!c             would yield a saving in time (TO DO if the selection is
!c             efficient)
!c
              do j=jcour+1,m
                  jol(j-1)=jol(j)
                  do i=1,n
                      ybar(i,j-1)=ybar(i,j)
                      sbar(i,j-1)=sbar(i,j)
                  enddo
              enddo
              jol(m)=iiter
              jcour=m
          end if
          ol(jcour)=olf
!c         print *,"jol=",(jol(i),i=1,m)
!c         write (6,'(a,i0,a,i0)') "nys=", nys, "  jcour=", jcour
!c
!c     --- selection by the Rayleigh quotient
!c
      elseif (select.eq.2) then
!c
!c>>>>>>>> the following should be changed in the futur so that a
!c         matrix build in a previous CG run could be updated (TO DO)
!c
!c         if (iiter.eq.1) nys=0
          if (nys.le.0) then
              nys=0
              jmin=1
              jmax=0
              do i=1,m
                  jol(i)=0
              enddo
          end if
!c
!c         --- it remains to set nys, jmax, jcour, ol(), and jol()
!c
          if (m.eq.1) then
              nys=1
              jmax=1
              jcour=1
              ol(1)=olf
              jol(1)=1
          elseif (nys.lt.m) then
!c
!c             --- select the current pair, since there is enough space
!c                 to save it
!c
              jmax=jmax+1
              jcour=jmax
              ol(jcour)=olf
!c
!c             --- it remains to set jol(); for this, sort the saved
!c                 pairs according to their OL factor
!c
              if (nys.eq.0) then
                  jol(1)=1
              else
!c
!c                 --- find the smallest index i2 such that
!c                     ol(jol(i2)) >= olf, this is the position of the
!c                     current OLF
!c
                  i2=1
  101             continue
                  if (ol(jol(i2)).ge.olf) goto 102
                      i2=i2+1
                      if (i2.le.nys) goto 101
  102             continue
!c
!c                 --- shift jol(i), for i>=i2
!c
                  if (i2.le.nys) then
                      do i=nys,i2,-1
                          jol(i+1)=jol(i)
                      enddo
                  end if
                  jol(i2)=jcour
              end if
              nys=nys+1
!c
          else
!c
!c             --- here, if the current (y,s) is selected, an older one
!c                 must be discarded;
!c
!c                 here, nys>=m (updated at the end of the ELSE), jmax=m
!c                 (no longer udated); it remains to set jcour (inside
!c                 the long if-then-else below and ol(jcour)=olf (at the
!c                 end of the ELSE);
!c
!c                 start by finding the smallest index i2 such that
!c                 ol(jol(i2)) >= olf, this is the position of the
!c                 current OLF
!c
              i2=1
  103         continue
              if (ol(jol(i2)).ge.olf) goto 104
                  i2=i2+1
                  if (i2.le.m) goto 103
  104         continue
              if (i2.eq.1) then
!c
!c                 --- olf is smaller (<=) than all the other OL factors
!c                     ==> save the (y,s) pair
!c
                  jcour=jol(1)
              elseif (i2.gt.m) then
!c
!c                 --- olf is greater (>) than all the other OL factors
!c                     ==> save the (y,s) pair
!c
                  jcour=jol(m)
              elseif (m.eq.2) then
!c
!c                 --- reject the (y,s) pair, since when m=2, we want to
!c                     save the pairs with the extreme OL factors
!c
                  mode=-1
                  return
              elseif (i2.eq.2) then
!c
!c                 --- check whether to replace the pair jol(2), by
!c                     looking at jol(1), jol(2), and jol(3)
!c
                  olmed=sqrt(ol(jol(1))*ol(jol(3)))
                  oldq=ol(jol(2))/olmed
                  if (oldq.lt.1.d0) oldq=1.d0/oldq
                  newq=olf/olmed
                  if (newq.lt.1.d0) newq=1.d0/newq
                  if (newq.ge.oldq) then
                      mode=-1
                      return
                  end if
                  jcour=jol(2)
              elseif (i2.eq.m) then
!c
!c                 --- check whether to replace the pair jol(m-1), by
!c                     looking at jol(m-2), jol(m-1), and jol(m)
!c
                  olmed=sqrt(ol(jol(m-2))*ol(jol(m)))
                  oldq=ol(jol(m-1))/olmed
                  if (oldq.lt.1.d0) oldq=1.d0/oldq
                  newq=olf/olmed
                  if (newq.lt.1.d0) newq=1.d0/newq
                  if (newq.ge.oldq) then
                      mode=-1
                      return
                  end if
                  jcour=jol(m-1)
              else
!c
!c                 --- here m>=4;
!c                     check whether to replace the pair jol(i2-1) or
!c                     jol(i2), by looking at jol(i2-2), jol(i2-1),
!c                     jol(i2), and jol(i2+1)
!c
                  ol0=log(ol(jol(i2-2)))
                  ol1=log(ol(jol(i2-1)))
                  ol2=log(ol(jol(i2)))
                  ol3=log(ol(jol(i2+1)))
                  olx=log(olf)
                  olvar =max(ol1-ol0,ol2-ol1,ol3-ol2)
                  olvar1=max(olx-ol0,ol2-olx,ol3-ol2)
                  olvar2=max(ol1-ol0,olx-ol1,ol3-olx)
                  if (olvar.le.min(olvar1,olvar2)) then
                      mode=-1
                      return
                  elseif (olvar1.lt.olvar2) then
                      jcour=jol(i2-1)
                  else
                      jcour=jol(i2)
                  end if
              end if
!c
!c             --- save the pair
!c
              nys=nys+1
              ol(jcour)=olf
          end if
!c         write (6,'(a,1pe13.6)') "olf = ", olf
!c         write (6,'(a)') " ol(jol(.))   jol(.)     ol(.)"
!c         write (6,'(1pe13.6,2x,i6,2x,1pe13.6)')
!c    &      (ol(jol(i)),jol(i),ol(i),i=1,min(nys,m))
!c         write (6,'(a,i0,a,i0)') "nys=", nys, "  jcour=", jcour
      else
!c
!c     --- unknown value of select
!c
          mode=3
          if (plev.ge.1) write (iom,903) select
  903     format (/" dysave-ERROR: unknown selection procedure,"," select = ",i5)
          return
      end if
!c
!c --- store ybar et sbar
!c
      r=sqrt(1.d0/ys)
      do i=1,n
          ybar(i,jcour)=r*y(i)
          sbar(i,jcour)=r*s(i)
      enddo
!c
!c --- save the OL factor
!c
!c     --- the simplest strategy is to take the OL factor at the
!c         first CG iteration (the function is supposed to be
!c         quadratic)
!c
      if (nys.eq.1) size=olf
!c
!c     --- another strategy is to take the geometric means of the
!c         previous OL factors (taking the OL factor, the closest
!c         in the logarithmic scale to this geometric means give
!c         approximately the same results)
!c
      olmed=1.d0
      do j=1,min(nys,m)
          olmed=olmed*ol(j)
      enddo
      olmed=olmed**(1.d0/dble(min(nys,m)))
      size=olmed
!c      print *,"ol=",(ol(j),j=1,min(nys,m))
!c      print *,"olmed=",olmed
!c
 1000 return
      end
!c
!c--------0---------0---------0---------0---------0---------0---------0--
!c
      subroutine deuclid (n,x,y,sp,izs,rzs,dzs)
      USE constants,    ONLY : dp
!c
      integer n,izs(1)
      real rzs(1)
      real(kind=dp) ::  x(n),y(n),sp,dzs(1)
!c
!c --- local variables
!c
      integer i
!c
      sp=0.d0
      do 10 i=1,n
          sp=sp+x(i)*y(i)
   10 continue
!c
      return
      end
!c
!c--------0---------0---------0---------0---------0---------0---------0--
!c

!c----------------------------------------------------------------------c
!c                                                                      c
!c     Product Mv                                                       c
!c                                                                      c
!c----------------------------------------------------------------------c
!c
      subroutine mvquad (n,a,v,av,izs,rzs,dzs)

  USE constants,                ONLY : dp , sp
!c
      integer ln
      common /cquad/ln
!c
      integer n,izs(1)
      real(kind=sp) ::  rzs(1)
      real(kind=dp) :: a(ln,n),v(n),av(n),dzs(1)
!c
      integer i,j
      real(kind=dp) :: r
!c
!c --- compute the matrix-vector product Av
!c
      do i=1,n
        r=0.d0
        do j=1,n
          r=r+a(i,j)*v(j)
        enddo
        av(i)=r
      enddo
!c
      return
      end



